
######################################################
## Replication file for Aidt and Leon (2020)        ##
## "The role of structural factors and diffusion    ##
## in social unrest: evidence from the Swing riots" ##
######################################################


rm(list=ls())

.libPaths("C:/Users/k1466170/Documents/R libraries")
.libPaths()

# Update as necessary
setwd("C:/Users/k1466170/Documents/Great Reform riots/UK riots/Data/Parish Names/Hobsbawm/replication material - BJPolS")

library(spdep)
library(foreign)
library(stargazer)
library(sphet)
library(plm)
library(glmmML)
library(splm)
library(jtools)
library(coefplot)
library(tidyverse)
library(dotwhisker)
library(gapminder)



###############
## Main Text ##
###############


# Load .dta file saved in old version (Stata 12 or earlier)
riots <- read.dta("riotscrosssectionRJan2020.dta")

# Create a matrix of locations
loc<-cbind(riots$longitude,riots$latitude)

# Find all neighbours within a 10 km radius (Great Circle distances)
nearby<-dnearneigh(loc,0,10,longlat=TRUE)

# Create the weight matrix
listw <- nb2listw(nearby,glist=NULL,style="W",zero.policy=TRUE)

# Create an index
id1 <- seq(1,nrow(loc))

# This creates the distances
tmp<-distance(loc, region.id = id1, output = TRUE, type = "distance", measure = "euclidean", shape.name="shapefile",
              region.id.name="id1", firstline=TRUE, file.name="dist10W.GWT")
rm(tmp)


coldist <- read.gwt2dist(file="dist10W.GWT", region.id=id1, skip=1)
class(coldist)
summary(coldist)
gc()


# MAIN RESULTS



# Table 1, column 1: SHAC with structural factors
res_T1C1 <- stslshac(allriotsX ~ cereal + ln_rural + enclosure + ln_wealth + ln_males + ln_total_of_person + ln_trade + ln_professionals + ln_petitions + news_10 + ln_garr_min_dis + Dpolice10 + factor(county),
                 data=riots, listw, na.action=na.fail, distance=coldist, zero.policy=TRUE, type="Triangular", HAC=TRUE)
summary(res_T1C1)

# Table 1, column 2: Run Swing replication.do

# Table 1, column 3: SHAC with structural factors and structural factors of neighbours (abridged version in text, full version in appendix CHECK)
res_T1C3 <- stslshac(allriotsX ~ cereal + ln_rural + enclosure + ln_wealth + ln_males + ln_total_of_person + ln_trade + ln_professionals + ln_petitions + news_10 + ln_garr_min_dis + Dpolice10 + Ccereal_10 + Cln_rural_10 + Cenclosure_10 + Cln_wealth_10 + Cln_males_10 + Cln_total_of_person_10 + Cln_trade_10 + Cln_professionals_10 + Cln_petitions_10 + Cnews10_10 + Cln_garr_min_dis_10 + CDpolice10_10 + factor(county),
                  data=riots, listw, na.action=na.fail, distance=coldist, zero.policy=TRUE, type="Triangular", HAC=TRUE)
summary(res_T1C3)


# Figure 2: Generated in file Figure 2.R



# Table 2

# Create the interactions
riots$int1 <- riots$cereal * riots$allriots_nearX
riots$int2 <- riots$ln_rural * riots$allriots_nearX
riots$int3 <- riots$enclosure * riots$allriots_nearX
riots$int4 <- riots$ln_wealth * riots$allriots_nearX
riots$int5 <- riots$ln_males * riots$allriots_nearX
riots$int6 <- riots$ln_total_of_person * riots$allriots_nearX
riots$int7 <- riots$ln_trade * riots$allriots_nearX
riots$int8 <- riots$ln_professionals * riots$allriots_nearX
riots$int9 <- riots$ln_petitions * riots$allriots_nearX
riots$int10 <- riots$news_10 * riots$allriots_nearX
riots$int11 <- riots$ln_garr_min_dis * riots$allriots_nearX
riots$int12 <- riots$Dpolice10 * riots$allriots_nearX


res_T2C1 <- stslshac(allriotsX ~ cereal + int1 + ln_rural + int2 + enclosure +int3 + ln_wealth + int4 + ln_males + int5 + ln_total_of_person + int6 + ln_trade + int7 + ln_professionals + int8 + ln_petitions + int9 + news_10 + int10 + ln_garr_min_dis + int11 + Dpolice10 + int12 + factor(county),
                  data=riots, listw, na.action=na.fail, distance=coldist, zero.policy=TRUE, type="Triangular", HAC=TRUE)
summary(res_T2C1, infer = c(TRUE, TRUE), level = .95)




##############
## APPENDIX ##
##############


# Figure A1: Distribution of distances to newspapers.
# Generated in newspaper distribution.do

# Table A2: Correlations
# Generated in Swing replication.do

# Table A3: Summary statistics
# Generated in Swing replication.do

# Table A4: The effect of structural factors of neighbors
# Generated above in Table 2

# Table A5: Core regression but with only small and only large riot

# Large riots only:
res_TA5C1 <- stslshac(allriots_largeX ~ cereal + ln_rural + enclosure + ln_wealth + ln_males + ln_total_of_person + ln_trade + ln_professionals + ln_petitions + news_10 + ln_garr_min_dis + Dpolice10 + factor(county),
                  data=riots, listw, na.action=na.fail, distance=coldist, zero.policy=TRUE, type="Triangular", HAC=TRUE)
summary(res_TA5C1)

# Small riots only:
res_TA5C2 <- stslshac(allriots_smallX ~ cereal + ln_rural + enclosure + ln_wealth + ln_males + ln_total_of_person + ln_trade + ln_professionals + ln_petitions + news_10 + ln_garr_min_dis + Dpolice10 + factor(county),
                       data=riots, listw, na.action=na.fail, distance=coldist, zero.policy=TRUE, type="Triangular", HAC=TRUE)
summary(res_TA5C2)


# Table A6: Robustness checks

# Column 1: 20 km neighborhood
# This requires a new weight matrix

# Create a matrix of locations
loc20<-cbind(riots$longitude,riots$latitude)
# Find all neighbours within a 20 km radius (Great Circle distances)
nearby20<-dnearneigh(loc20,0,20,longlat=TRUE)

# Create the weight matrix
listw20 <- nb2listw(nearby20,glist=NULL,style="W",zero.policy=TRUE)

# Create an index
id20 <- seq(1,nrow(loc20))

# This creates the distances
tmp20<-distance(loc20, region.id = id20, output = TRUE, type = "distance", measure = "euclidean", shape.name="shapefile",
                region.id.name="id20", firstline=TRUE, file.name="dist20W.GWT")
rm(tmp20)

coldist20 <- read.gwt2dist(file="dist20W.GWT", region.id=id20, skip=1)
class(coldist20)
summary(coldist20)
gc()

res_A6C1 <- stslshac(allriotsX ~ cereal + ln_rural + enclosure + ln_wealth + ln_males+ ln_total_of_person + ln_trade + ln_professionals + ln_petitions + news_10 + ln_garr_min_dis + Dpolice10 + factor(county),
                 data=riots, listw20, na.action=na.fail, distance=coldist20, zero.policy=TRUE, type="Triangular", HAC=TRUE)
summary(res_A6C1)



# COLUMN 2: 30 km neighbourhood
# This requires a new weight matrix

# Create a matrix of locations
loc30<-cbind(riots$longitude,riots$latitude)
# Find all neighbours within a 30 km radius (Great Circle distances)
nearby30<-dnearneigh(loc30,0,30,longlat=TRUE)

# Create the weight matrix
listw30 <- nb2listw(nearby30,glist=NULL,style="W",zero.policy=TRUE)

# Create an index
id30 <- seq(1,nrow(loc30))

# This creates the distances
tmp30<-distance(loc30, region.id = id30, output = TRUE, type = "distance", measure = "euclidean", shape.name="shapefile",
                region.id.name="id30", firstline=TRUE, file.name="dist30W.GWT")
rm(tmp30)

coldist30 <- read.gwt2dist(file="dist30W.GWT", region.id=id30, skip=1)
class(coldist30)
summary(coldist30)
gc()

res_A4C2 <- stslshac(allriotsX ~ cereal + ln_rural + enclosure + ln_wealth + ln_males + ln_total_of_person + ln_trade + ln_professionals + ln_petitions + news_10 + ln_garr_min_dis + Dpolice10 + factor(county),
                  data=riots, listw30, na.action=na.fail, distance=coldist30, zero.policy=TRUE, type="Triangular", HAC=TRUE)
summary(res_A4C2)



# COLUMN 3: Without Kent
# This requires a new weight matrix

riotsNK <- subset (riots, DKent != 1)

# Create a matrix of locations
locNK<-cbind(riotsNK$longitude,riotsNK$latitude)

# Find all neighbours within a 10 km radius (Great Circle distances)
nearbyNK<-dnearneigh(locNK,0,10,longlat=TRUE)

# Create the weight matrix
#W<-nb2mat(nearby,glist=NULL,style="B",zero.policy=TRUE)
listwNK <- nb2listw(nearbyNK,glist=NULL,style="W",zero.policy=TRUE)

# Create an index
idNK <- seq(1,nrow(locNK))

# This creates the distances
tmpNK<-distance(locNK, region.id = idNK, output = TRUE, type = "distance", measure = "euclidean", shape.name="shapefile",
                region.id.name="idNK", firstline=TRUE, file.name="distNKW.GWT")
rm(tmpNK)


coldistNK <- read.gwt2dist(file="distNKW.GWT", region.id=idNK, skip=1)
class(coldistNK)
summary(coldistNK)
gc()

res_A4C3 <- stslshac(allriotsX ~cereal + ln_rural + enclosure + ln_wealth + ln_males + ln_total_of_person + ln_trade + ln_professionals + ln_petitions + news_10 + ln_garr_min_dis + Dpolice10 + factor(county),
                 data=riotsNK, listwNK, na.action=na.fail, distance=coldistNK, zero.policy=TRUE, type="Triangular", HAC=TRUE)
summary(res_A4C3)


# COLUMN 4: Without London
# This requires a new weight matrix

riotsNL <- subset (riots, DLondon20 != 1)

# Create a matrix of locations
locNL<-cbind(riotsNL$longitude,riotsNL$latitude)

# Find all neighbours within a 10 km radius (Great Circle distances)
nearbyNL<-dnearneigh(locNL,0,10,longlat=TRUE)

# Create the weight matrix
listwNL <- nb2listw(nearbyNL,glist=NULL,style="W",zero.policy=TRUE)

# Create an index
idNL <- seq(1,nrow(locNL))

# This creates the distances
tmpNL<-distance(locNL, region.id = idNL, output = TRUE, type = "distance", measure = "euclidean", shape.name="shapefile",
                region.id.name="idNL", firstline=TRUE, file.name="distNLW.GWT")
rm(tmpNL)

coldistNL <- read.gwt2dist(file="distNLW.GWT", region.id=idNL, skip=1)
class(coldistNL)
summary(coldistNL)
gc()

res_A4C4 <- stslshac(allriotsX ~ cereal + ln_rural + enclosure + ln_wealth + ln_males + ln_total_of_person + ln_trade + ln_professionals + ln_petitions + news_10 + ln_garr_min_dis + Dpolice10 + factor(county),
                 data=riotsNL, listwNL, na.action=na.fail, distance=coldistNL, zero.policy=TRUE, type="Triangular", HAC=TRUE)
summary(res_A4C4)

